/*
 * SystemState.cpp
 *
 *  Created on: 22 Aug 2011
 *      Author: Allan
 */

#include "SystemState.h"

// C++ includes
#include <cmath>

// GeoReact includes
#include "ElectrolyteSolution.h"
#include "Utilities/Algorithms.h"

SystemState::SystemState()
{}

SystemState::SystemState(const System& system) :
System(system)
{
}

SystemState::SystemState(const System& system, double T, double P, const VectorXd& n) :
System(system)
{
	SetState(T, P, n);
}

void SystemState::SetState(double T, double P, const VectorXd& n)
{
	this->T = T;
	this->P = P;
	this->n = n;
	
	a = multiphase.Activities(T, P, n);
	g = multiphase.ActivityCoefficients(T, P, n);
}

const double SystemState::Get(const string& quantity) const
{
	if(quantity == "IonicStrength" || quantity == "IonicStrength(Effective)")
	{
		const vector<string> aqueousSpecies = GetPhase("Aqueous").GetSpecies();
		
		ElectrolyteSolution sol(aqueousSpecies);
		
		return sol.EffectiveIonicStrength(GetPhaseComposition("Aqueous", n));
	}
	
	if(quantity == "IonicStrength(Stoichiometric)")
	{
		const vector<string> aqueousSpecies = GetPhase("Aqueous").GetSpecies();
		
		ElectrolyteSolution sol(aqueousSpecies);
		
		return sol.StoichiometricIonicStrength(GetPhaseComposition("Aqueous", n));
	}
	
	if(quantity == "pH")
	{
		const SpeciesMap& speciesMap = multiphase.GetSpeciesMap();
		
		double iH = speciesMap["H[+]"];
		
		return -std::log10(a[iH]);
	}
	
	cerr << "Error: The quantity type " << quantity << " is not valid." << endl; return 0.0;
}

const double SystemState::Get(const string& quantity, const string& species) const
{
	const SpeciesMap& speciesMap = multiphase.GetSpeciesMap();
	
	const Index ispecies = speciesMap[species];
	
	if(quantity == "Moles" || quantity == "Concentration")
		return n[ispecies];
	
	if(quantity == "Activity")
		return a[ispecies];
	
	if(quantity == "ActivityCoefficient")
		return g[ispecies];
	
	if(quantity == "PartialPressure")
	{
		const VectorXd ng = GetPhaseComposition("Gaseous", n);
		
		return P * (n[ispecies] / ng.sum());
	}
	
	if(quantity == "SI" || quantity == "SaturationIndex")
	{
		const Reaction& reaction = reactionMap.find(species)->second;
		
		return std::log10(reaction.ReactionQuotient(a) / reaction.GetEquilibriumConstant()(T, P));
	}
	
	if(quantity == "K" || quantity == "EquilibriumConstant")
	{
		const Reaction& reaction = reactionMap.find(species)->second;
		
		return reaction.GetEquilibriumConstant()(T, P);
	}
	
	if(quantity == "MassRate")
	{
		const Reaction& reaction = reactionMap.find(species)->second;

		return reaction.MassRate(T, P, n, a);
	}
	
	cerr << "Error: The quantity type " << quantity << " is not valid." << endl; return 0.0;
}

void SystemState::Print() const
{
	System::Print();
	
	cout << "Temperature: " << T << " Celsius" << endl;
	cout << "Pressure: " << P << " bar" << endl << endl;
	
	multiphase.PrintState(T, P, n);
	
	cout << endl << endl;
	
	eReactions.PrintState(T, P, n, a);
	
	cout << endl << endl;
	
	kReactions.PrintState(T, P, n, a);
}
